# Estimate our main model & define the weight matrix within Province.

# remove everything
rm(list=ls())

# =================================================================================
# ---------------------------------------------------------------------
# Drop outliers in dataset "Data" according to the "v"th column
# at each percentage "p" at bottom and top. The default value of
# p is 0.005. v is a vector of numbers.
# ---------------------------------------------------------------------

outlier <- function(Data, v, p=0.005) {
  out <- c()
  for(i in v) {
    out <- rbind(out, quantile(Data[,i], probs=c(p, 1-p)))
  }
  j <- 1
  for(i in v) {
    Data <- Data[Data[,i]>out[j,1] & Data[,i]<out[j,2], ]; print(dim(Data))
    j <- j+1
  }
  return(Data)
}

# =================================================================================

# Load packages
library(foreign)
library(gmm)

# read data 
Data <- read.dta('Ind39.dta')

Data <- outlier(Data, which(names(Data)=="sm"), p=0.005)

# Delete obs according to variable "Foreign", "Province", and "sm"
Data <- subset(Data, F1>=0 & F1<=1 & !is.na(Province) & sm<=1)
# summary(Data)

# =============================================
# Step One
h.lnbe <- mean(log(Data$sm)); h.eta <- - log(Data$sm) + h.lnbe; h.e <- mean(exp(h.eta))
h.betam <- exp(h.lnbe)/h.e

# =============================================
# Step Two
# Generate weights & variables related to weights
for(i in 1:nrow(Data)) {
  s1.t <- (Data$Year==Data$Year[i])*(Data$Province==Data$Province[i]); s1.t[i] <- 0
  s1.b <- sum(s1.t)
  s1 <- if(s1.b>0) s1.t/s1.b else rep(0, length(s1.t));
  Data$wk1[i] <- sum(s1*Data$k1); Data$wm1[i] <- sum(s1*Data$m1, na.rm = TRUE); Data$wl1[i] <- sum(s1*Data$l1, na.rm = TRUE)
}

# Generate year dummies
TD <- model.matrix(~as.factor(Data$Year)-1); TD1 <- model.matrix(~as.factor(Data$Year-1)-1)[,-1]

# Define functions: "wphi" and "phi"
wphi <- function(beta) {
  (1-h.betam)*Data$wm1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$wk1 - beta[2]*Data$wl1 - TD1%*%beta[3:10] 
}

phi <- function(beta) {
  (1-h.betam)*Data$m1 - h.lnbe - log(Data$ppi1/Data$ppii1) - beta[1]*Data$k1 - beta[2]*Data$l1 - TD1%*%beta[3:10]
}

# -----------------------------------------------
# Optimization

objfn <- function(beta){
  y.star <- Data$y - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
  p1 <- as.vector(phi(beta)); p2 <- as.vector(wphi(beta)); 
  var.poly <- poly(cbind(p1, Data$F1, p2), degree = 2, raw = TRUE)
  mod1 <- lm(as.vector(y.star)~var.poly)
  OBJ <- sum(mod1$residuals^2)
  return(OBJ)
}

set.seed(1)
initial <- c(0.05, 0.1, 0.05, 0,0,0,0,0,0,0,0,0)
target <- optim(initial, objfn, method = "Nelder-Mead")

# results - elas
h.betak <- target$par[1]; h.betal <- target$par[2]

G.coef <- function(beta){
  y.star <- Data$y - h.betam*Data$m - beta[1]*Data$k - beta[2]*Data$l - TD%*%beta[3:11]
  p1 <- as.vector(phi(beta)); p2 <- as.vector(wphi(beta)); 
  var.poly <- poly(cbind(p1, Data$F1, p2), degree = 2, raw = TRUE)
  mod1 <- lm(as.vector(y.star)~var.poly)
  return(mod1$coef)
}
# results - partial effects
coef <- G.coef(target$par)[2:10]
w1 <- as.vector(phi(target$par)); sw1 <- as.vector(wphi(target$par))

p.w  <- coef[1] + 2*coef[2]*w1 + coef[4]*Data$F1 + coef[7]*sw1
p.F  <- coef[3] + coef[4]*w1 + 2*coef[5]*Data$F1 + coef[8]*sw1
p.sw <- coef[6] + coef[7]*w1 + coef[8]*Data$F1 + 2*coef[9]*sw1

plot(density(p.F)); r2 <- summary(p.F)
plot(density(p.sw)); r3 <- summary(p.sw)

# The estimated w and w1
Data$w <- Data$y - h.betak*Data$k - h.betal*Data$l - h.betam*Data$m - TD%*%target$par[3:11] - h.eta
Data$w1 <- phi(target$par); Data$sw1 <- wphi(target$par); Data$eta <- h.eta

# Report the basic results
r00 <- c("","Min", "Q1","Median","Mean","Q3","Max")
r2 <- c("G", r2)
r3 <- c("sw", r3)
Table.t <- rbind(r00, r3, r2); rownames(Table.t) <- c(); colnames(Table.t) <- c()
print(Table.t, quote = FALSE)


# save data
Data$pF <- p.F; Data$pw <- p.w; Data$psw <- p.sw
Data39 <- Data[ , c("LPCode","Firm","Year","Province","F1","F2","wk1","wm1","wl1","w","w1","sw1","eta","pF","pw","psw")]
write.dta(Data39, file = "Ind39_s_Dt_Est.dta")
